MCMC convergence diagnostics
Adam Howes
Abstract
Background In the report naomi-simple_fit with parameter tmbstan = TRUE, we used the NUTS algorithm to perform MCMC inference for the simplified Naomi model.
Task Here we assess whether or not the results of the MCMC are suitable using a range of diagnostic tools.
We start by obtaining results from the latest version of naomi-simple_fit with tmbstan = TRUE.
out <- readRDS("depends/out.rds")
mcmc <- out$mcmcThis MCMC took 18.9085492 minutes to run
bayesplot::color_scheme_set("viridisA")
ggplot2::theme_set(theme_minimal())1 \(\hat R\)
We are looking for values of \(\hat R\) less than 1.1 here.
rhats <- bayesplot::rhat(mcmc)
bayesplot::mcmc_rhat(rhats)
(big_rhats <- rhats[rhats > 1.1])## beta_rho[1] beta_rho[2] beta_alpha[1] beta_alpha[2] logit_phi_rho_x log_sigma_rho_x logit_phi_rho_xs log_sigma_rho_xs logit_phi_rho_a log_sigma_rho_a logit_phi_rho_as log_sigma_rho_as u_rho_x[4]
## 1.436462 2.050365 1.273354 1.462126 1.248576 1.150130 1.494549 1.362138 1.126501 1.151890 1.256854 1.325939 1.171248
## us_rho_x[4] us_rho_x[6] u_rho_xs[1] u_rho_xs[2] u_rho_xs[4] u_rho_xs[11] u_rho_xs[12] u_rho_xs[22] u_rho_xs[27] u_rho_xs[30] us_rho_xs[1] us_rho_xs[2] us_rho_xs[3]
## 1.139724 1.159936 1.108818 1.104644 1.118410 1.113635 1.110012 1.174033 1.201293 1.139680 1.127415 1.138374 1.174318
## us_rho_xs[4] us_rho_xs[5] us_rho_xs[6] us_rho_xs[8] us_rho_xs[9] us_rho_xs[10] us_rho_xs[11] us_rho_xs[12] us_rho_xs[13] us_rho_xs[14] us_rho_xs[15] us_rho_xs[22] us_rho_xs[23]
## 1.133138 1.156225 1.183878 1.133482 1.202853 1.187893 1.199672 1.198588 1.185517 1.125922 1.194364 1.185790 1.110124
## us_rho_xs[24] us_rho_xs[25] us_rho_xs[26] us_rho_xs[27] us_rho_xs[28] us_rho_xs[29] us_rho_xs[30] us_rho_xs[31] us_rho_xs[32] u_rho_a[1] u_rho_a[2] u_rho_a[3] u_rho_a[4]
## 1.238823 1.134849 1.259984 1.296719 1.176014 1.247366 1.324837 1.255172 1.188287 1.429221 1.429893 1.435065 1.430681
## u_rho_a[5] u_rho_a[6] u_rho_a[7] u_rho_a[8] u_rho_a[9] u_rho_a[10] u_rho_as[1] u_rho_as[2] u_rho_as[3] u_rho_as[4] u_rho_as[5] u_rho_as[6] u_rho_as[7]
## 1.432134 1.432472 1.432241 1.436418 1.434241 1.434673 2.017298 2.033510 2.041822 2.038281 2.054323 2.044709 2.034816
## u_rho_as[8] u_rho_as[9] u_rho_as[10] logit_phi_alpha_x log_sigma_alpha_x logit_phi_alpha_xs log_sigma_alpha_xs u_alpha_x[1] u_alpha_x[2] u_alpha_x[3] u_alpha_x[4] u_alpha_x[8] u_alpha_x[9]
## 2.029767 2.036252 1.996995 1.101527 1.410485 1.146635 1.268746 1.263213 1.166815 1.240560 1.177752 1.193106 1.171466
## u_alpha_x[10] u_alpha_x[13] u_alpha_x[14] u_alpha_x[21] u_alpha_x[22] u_alpha_x[23] u_alpha_x[24] u_alpha_x[25] u_alpha_x[26] u_alpha_x[27] u_alpha_x[28] u_alpha_x[30] us_alpha_x[1]
## 1.187447 1.273750 1.117343 1.432030 1.195067 1.190772 1.272183 1.363360 1.142784 1.178192 1.109855 1.130643 1.168862
## us_alpha_x[2] us_alpha_x[3] us_alpha_x[4] us_alpha_x[5] us_alpha_x[6] us_alpha_x[8] us_alpha_x[12] us_alpha_x[15] us_alpha_x[21] us_alpha_x[22] us_alpha_x[24] us_alpha_x[25] us_alpha_x[29]
## 1.133704 1.197043 1.264392 1.112401 1.198646 1.100455 1.118491 1.118520 1.189470 1.245977 1.164936 1.132086 1.162871
## us_alpha_x[30] u_alpha_xs[1] u_alpha_xs[8] u_alpha_xs[11] u_alpha_xs[13] u_alpha_xs[16] u_alpha_xs[17] u_alpha_xs[18] u_alpha_xs[19] u_alpha_xs[21] u_alpha_xs[22] u_alpha_xs[23] u_alpha_xs[24]
## 1.136288 1.115010 1.139084 1.103883 1.162435 1.145249 1.133684 1.299266 1.209079 1.405892 1.201692 1.144698 1.275732
## u_alpha_xs[25] u_alpha_xs[28] u_alpha_xs[29] u_alpha_xs[32] us_alpha_xs[1] us_alpha_xs[2] us_alpha_xs[3] us_alpha_xs[4] us_alpha_xs[5] us_alpha_xs[6] us_alpha_xs[8] us_alpha_xs[9] us_alpha_xs[10]
## 1.257591 1.187026 1.390206 1.111991 1.158874 1.175500 1.192279 1.138272 1.181611 1.150340 1.217743 1.187316 1.192677
## us_alpha_xs[11] us_alpha_xs[12] us_alpha_xs[13] us_alpha_xs[15] us_alpha_xs[17] us_alpha_xs[18] us_alpha_xs[19] us_alpha_xs[20] us_alpha_xs[21] us_alpha_xs[22] us_alpha_xs[23] us_alpha_xs[24] us_alpha_xs[25]
## 1.180764 1.122971 1.176139 1.124448 1.145104 1.138565 1.131682 1.104125 1.154424 1.203935 1.146553 1.133402 1.145399
## us_alpha_xs[26] u_alpha_a[1] u_alpha_a[2] u_alpha_a[3] u_alpha_a[4] u_alpha_a[5] u_alpha_a[6] u_alpha_a[7] u_alpha_a[8] u_alpha_a[9] u_alpha_a[10] u_alpha_a[11] u_alpha_a[12]
## 1.115193 1.237430 1.260247 1.265380 1.224632 1.234161 1.248809 1.257960 1.259736 1.262103 1.253853 1.255800 1.234374
## u_alpha_a[13] u_alpha_as[1] u_alpha_as[2] u_alpha_as[3] u_alpha_as[4] u_alpha_as[5] u_alpha_as[6] u_alpha_as[7] u_alpha_as[8] u_alpha_as[9] u_alpha_as[10] u_alpha_xa[13] u_alpha_xa[21]
## 1.218461 1.405071 1.431038 1.426739 1.435003 1.462654 1.469958 1.488702 1.483664 1.495276 1.430269 1.107958 1.335780
## u_alpha_xa[24] u_alpha_xa[25] u_alpha_xa[29] log_sigma_lambda_x ui_lambda_x[1] ui_lambda_x[10] ui_lambda_x[11] ui_lambda_x[13] ui_lambda_x[14] ui_lambda_x[22] ui_lambda_x[23] ui_lambda_x[32] log_sigma_ancalpha_x
## 1.259357 1.108961 1.204265 2.309844 1.212403 1.427155 1.258550 1.187641 1.158190 1.139450 1.131882 1.146537 1.470755
## ui_anc_rho_x[15] ui_anc_rho_x[25] ui_anc_rho_x[29] ui_anc_alpha_x[16] ui_anc_alpha_x[18] ui_anc_alpha_x[19] ui_anc_alpha_x[24] ui_anc_alpha_x[25] ui_anc_alpha_x[29] ui_anc_alpha_x[32] log_or_gamma[3] log_or_gamma[5] log_or_gamma[8]
## 1.101091 1.172969 1.114979 1.195326 1.303179 1.144301 1.125907 1.151146 1.269189 1.129052 1.124149 1.113333 1.101212
## log_or_gamma[9] log_or_gamma[10] log_or_gamma[12] log_or_gamma[18] log_or_gamma[20] log_or_gamma[22] log_or_gamma[31] lp__
## 1.123776 1.127447 1.141006 1.102396 1.114067 1.128760 1.113397 1.567209length(big_rhats) / length(rhats)## [1] 0.41260162 ESS ratio
Reasonable to be worried about values less than 0.1 here.
ratios <- bayesplot::neff_ratio(mcmc)
bayesplot::mcmc_neff(ratios)
3 Autocorrelation
How much autocorrelation is there in the chains?
bayesplot::mcmc_acf(mcmc, pars = vars(starts_with("beta")))
4 Univariate traceplots
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("beta")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("logit")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_sigma")))
4.1 Prevalence model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_a[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_as[")))
4.2 ART model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_a[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_as[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xa[")))
4.3 Other
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_lambda_x[")))
4.4 ANC model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_or_gamma["))) #' N.B. these are from the ANC attendance model
5 Pairs plots
There is a prior suspicion (from Jeff, Tim, Rachel) that the ART attendance model is unidentifiable.
Let’s have a look at the pairs plot for neighbouring districts and the log_or_gamma parameter.
nb <- area_merged %>%
filter(area_level == max(area_level)) %>%
bsae::sf_to_nb()
neighbours_pairs_plot <- function(par, i) {
neighbour_pars <- paste0(par, "[", c(i, nb[[i]]), "]")
bayesplot::mcmc_pairs(mcmc, pars = neighbour_pars, diag_fun = "hist", off_diag_fun = "hex")
}
# area_merged %>%
# filter(area_level == max(area_level)) %>%
# print(n = Inf)Here are Nkhata Bay and neighbours:
neighbours_pairs_plot("log_or_gamma", 5) 
And here are Blantyre and neighbours:
neighbours_pairs_plot("log_or_gamma", 26)
6 NUTS specific assessment
np <- bayesplot::nuts_params(mcmc)Are there any divergent transitions?
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_divergent = sum(Value))bayesplot::mcmc_nuts_divergence(np, bayesplot::log_posterior(mcmc))
We can also use energy plots (Betancourt 2017): ideally these two histograms would be the same When the histograms are quite different, it may suggest the chains are not fully exploring the tails of the target distribution.
bayesplot::mcmc_nuts_energy(np)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.